import matplotlib.pyplot as plt
import cobra
from cobra.io import validate_sbml_model
import importlib
from xml.etree import ElementTree
import utils.Model_correction as mc
import sys
import utils.model_maj as mj
import utils.viz_utils as vu
#import cplex
import cbmpy
import pandas as pd
import plotly.express as px
pyparsing import INFO: No xlwt module available, Excel spreadsheet creation disabled CBGLPK based on swiglpk: not all methods implimented yet! 5.0 ***** Using CPLEX ***** doFBAMinSum not available with GLPK INFO: No xlrd module available, Excel spreadsheet reading disabled *********************************************************************** * Welcome to CBMPy (0.8.4) - PySCeS Constraint Based Modelling * * http://cbmpy.sourceforge.net * * Copyright(C) Brett G. Olivier 2014 - 2020 * * Systems Biology Lab, Vrije Universiteit Amsterdam * * Amsterdam, The Netherlands * * CBMPy is developed as part of the BeBasic MetaToolKit Project * * Distributed under the GNU GPL v 3.0 licence, see * * LICENCE (supplied with this release) for details * ***********************************************************************
importlib.reload(mj)
importlib.reload(vu)
importlib.reload(mc)
<module 'utils.Model_correction' from '/home/jmuller/Documents/Stage_metabo/code/utils/Model_correction.py'>
HepG2, errors = validate_sbml_model("../models_storage/Hep-G2_v3.xml")
iHep, errors = validate_sbml_model("../models_storage/iHepatocytes2322_Cobra.xml")
HepG2 = mj.maj("iHepatocytes2322_Cobra.xml", "Hep-G2_v2.xml", "../models_storage/", bounds_check = False, genes_id_copy= False, alt_gene_ids= False, metab_id_check= False, bounds_value_check= False, subsystem_copy= True )
# Bidouillage pour ajouter des subsystems aux réactions d'échange. (Va ajouter "exchange reactions" normalement.)
for r in HepG2.reactions :
if "EX_" in r.id :
m = [m_id for m_id, _ in HepG2.reactions.get_by_id(r.id).metabolites.items()]
sub = [r.subsystem for r in m[0].reactions if "EX_" not in r.id]
r.subsystem = sub[0]
HepG2.reactions.EX_m01965x.bounds = (-1000.0,-1000.0) #Glucose exchange
HepG2.reactions.HMR_3883.bounds = (0.0,0.0) #Serine --> Pyruvate exchange
HepG2.reactions.EX_m02630x.bounds = (-1000.0,1000.0) #O2 exchange
HepG2.reactions.EX_m02819x.bounds = (0.0,1000.0) #Pyruvate exchange
HepG2.reactions.EX_m01910x.bounds = (0.0,0.0) #Galactose exchange
HepG2.reactions.EX_m01743x.bounds = (0.0, 1000.0) #D-Ribulose exchange
HepG2.reactions.EX_m01962x.bounds = (0.0,1000.0) #glucosamine exchange
HepG2.reactions.HMR_4316.bounds = (-1000.0,1000.0) # Glucose --> D-Glucitol
HepG2.reactions.EX_m02896x.bounds = (0.0,1000.0) # Serine intake
HepG2.reactions.EX_m01682x.bounds = (-1000.0,0.0) # Glucitol secretion blocked. Kept the intake just in case.
HepG2.reactions.EX_m01840x.bounds = (-1000.0,0.0) #Fructose exchange
#HepG2.reactions.HMR_4930.bounds = (-1000.0,1000.0) # Pyruvate transfer from cytoplasm to peroxysome
#HepG2.reactions.HMR_4281.bounds = (0.0,0.0) # Fermentation in peroxysome
#HepG2.reactions.EX_m02403x.bounds = (0.0,1000.0) #Lactate exchange
HepG2.reactions.EX_m02819x.bounds = (0.0,1000.0) #Pyruvate exchange
#HepG2.reactions.EX_m01840x.bounds = (0.0,0.0) #Fructose exchange
HepG2.reactions.HMR_4281.bounds = (-1000.0,0.0) # Pyruvate fermentation in peroxysome
HepG2.reactions.EX_m02630x.bounds = (0.0,1000.0) #O2 exchange
HepG2.reactions.EX_m02630x.bounds = (-1000.0,1000.0) #O2 exchange
HepG2.reactions.EX_m02819x.bounds = (0.0,1000.0) #Pyruvate exchange
HepG2.reactions.EX_m01910x.bounds = (0.0,0.0) #Galactose exchange
HepG2.reactions.EX_m01743x.bounds = (0.0, 1000.0) #D-Ribulose exchange
HepG2.reactions.EX_m01962x.bounds = (0.0,1000.0) #glucosamine exchange
#HepG2.reactions.EX_m01965x.bounds = (-1000,1000.0) #Glucose exchange
#HepG2.reactions.EX_m01286x.bounds = (0.0,0.0) #ADP-Glucose exchange
#HepG2.reactions.EX_m01840x.bounds = (-1000.0,0.0) #Fructose exchange
#HepG2.reactions.EX_m01743x.bounds = (-1000.0,1000.0) #D-Ribulose exchange
#HepG2.reactions.EX_m02453x.bounds = (0.0,0.0) #Mannose exchange
#HepG2.reactions.EX_m02470x.bounds = (0.0,0.0) #Methanol exchange
# --> After Methanol cut, Glucose uptake flux went from 0.0 to -190.0
# Without O2, glucose uptake went down to -23
HepG2.objective = "biomass_components"
sol = HepG2.optimize()
sol
| fluxes | reduced_costs | |
|---|---|---|
| EX_m00097x | -0.000000 | 0.000000e+00 |
| EX_m00157x | 952.717900 | 0.000000e+00 |
| EX_m00228x | 0.000000 | 0.000000e+00 |
| EX_m00242x | 0.000000 | 0.000000e+00 |
| EX_m00266x | 0.000000 | 0.000000e+00 |
| ... | ... | ... |
| HMR_9715 | 0.000000 | 0.000000e+00 |
| HMR_9730 | 0.000000 | 0.000000e+00 |
| HMR_9736 | 36.365246 | 0.000000e+00 |
| biomass_components | 27.150400 | 4.163336e-17 |
| artificial_biomass | 0.000000 | -2.000000e+00 |
5906 rows × 2 columns
new_iHep.objective = "biomass_components"
sol = new_iHep.optimize()
sol
[(m.name, m.id) for m in HepG2.reactions.biomass_components.metabolites]
importlib.reload(vu)
vu.print_exchanges(HepG2, filter = "all")
f = vu.run_parcours(HepG2.reactions.EX_m01965x)
print(len(f))
#m01682c glucitol
for reaction, flux in f.items() :
vu.print_reactions(HepG2.reactions.get_by_id(reaction), flux)
print(f"FLUX : {flux} --- ID : {HepG2.reactions.get_by_id(reaction).id} --- compartment : {HepG2.reactions.get_by_id(reaction).compartments} \n\n---\n\n")
vu.print_reactions(HepG2.reactions.biomass_components)
from cobra.flux_analysis import flux_variability_analysis
fva = flux_variability_analysis(HepG2)
indispensable = []
unisens = []
reversible = []
unknown = []
for reaction_id, flux_values in fva.iterrows() :
min_flux = flux_values["minimum"]
max_flux = flux_values["maximum"]
difference = abs(max_flux - min_flux)
if max_flux >0.0 and min_flux > 0.0 :
f = 1.0
elif max_flux > 0.0 and min_flux < 0.0 :
f = 0.0
elif max_flux < 0.0 and min_flux < 0.0 :
f = -1.0
else :
f = 0.0
if difference < 0.01 and min_flux + max_flux != 0.0 and abs(min_flux + max_flux) > 0.01 :
indispensable.append((reaction_id, flux_values))
elif min_flux < 0.0 and max_flux > 0.0 :
reversible.append((reaction_id, flux_values))
elif min_flux > 0.0 and max_flux > 0.0 or min_flux < 0.0 and max_flux < 0.0 :
if difference > 0.01 :
unisens.append((reaction_id, flux_values))
else :
unknown.append((reaction_id, flux_values))
print("\n--------INDISPENSABLE--------\n")
for reaction_id, flux_values in indispensable :
min_flux = flux_values["minimum"]
max_flux = flux_values["maximum"]
vu.print_reactions(HepG2.reactions.get_by_id(reaction_id), f)
print(f"\nMIN FLUX -- {min_flux} === {max_flux} -- MAX FLUX\n\n")
print("\n--------UNISENS--------\n")
for reaction_id, flux_values in unisens :
min_flux = flux_values["minimum"]
max_flux = flux_values["maximum"]
vu.print_reactions(HepG2.reactions.get_by_id(reaction_id), f)
print(f"\nMIN FLUX -- {min_flux} === {max_flux} -- MAX FLUX\n\n")
print("\n--------REVERSIBLE--------\n")
for reaction_id, flux_values in reversible :
min_flux = flux_values["minimum"]
max_flux = flux_values["maximum"]
vu.print_reactions(HepG2.reactions.get_by_id(reaction_id), f)
print(f"\nMIN FLUX -- {min_flux} === {max_flux} -- MAX FLUX\n\n")
if len(unknown) >0 :
print("\n--------UNKNOWN--------\n")
for reaction_id, flux_values in unknown :
min_flux = flux_values["minimum"]
max_flux = flux_values["maximum"]
vu.print_reactions(HepG2.reactions.get_by_id(reaction_id), f)
print(f"\nMIN FLUX -- {min_flux} === {max_flux} -- MAX FLUX\n\n")
df_C, df_R, df_S, df_M, df_P, df_X, df_L, df_G, df_N = vu.build_reaction_df(HepG2)
flux_filter = 0.0
fig = px.treemap(df_C.loc[(df_C["flux"] >= flux_filter)] , path=['subSystem', 'id'],
values='flux', color='subSystem',color_discrete_sequence=px.colors.qualitative.Light24_r)
fig.update_layout(title_text="Cytoplasm reactions", font_size=12)
fig.show(renderer='notebook')
fig = px.treemap(df_R.loc[df_R["flux"] >= flux_filter], path=['subSystem', 'id'],
values='flux', color='subSystem',color_discrete_sequence=px.colors.qualitative.Light24_r)
fig.update_layout(title_text="Endoplasmic Reticulum reactions", font_size=12)
fig.show(renderer='notebook')
fig = px.treemap(df_S.loc[df_S["flux"] >= flux_filter], path= ['subSystem', 'id'],
values='flux', color='subSystem',color_discrete_sequence=px.colors.qualitative.Light24_r)
fig.update_layout(title_text="Endomembraneous reactions", font_size=12)
fig.show(renderer='notebook')
fig = px.treemap(df_M.loc[df_M["flux"] >= flux_filter], path= [px.Constant("Mitochondria"),'subSystem', 'id'],
values='flux', color='subSystem',color_discrete_sequence=px.colors.qualitative.Light24_r)
fig.update_layout(title_text="Mitochondrial reactions", font_size=12)
fig.show(renderer='notebook')
fig = px.treemap(df_P.loc[df_P["flux"] >= flux_filter], path= [px.Constant("Peroxysome"),'subSystem', 'id'],
values='flux', color='subSystem',color_discrete_sequence=px.colors.qualitative.Light24_r)
fig.update_layout(title_text="Peroxysomal reactions", font_size=12)
fig.show(renderer='notebook')
fig = px.treemap(df_X.loc[df_X["flux"] >= flux_filter], path= ['subSystem', 'id'],
values='flux', color='subSystem',color_discrete_sequence=px.colors.qualitative.Light24_r)
fig.update_layout(title_text="Exchange reactions", font_size=12)
fig.show(renderer='notebook')